## To create and send PBS files to run natural waterbody simulations ###
# as for spread_model_Pilbara_barriers_PBS.R except we imagine that toads can move for 180 days

setwd("~/final_natural_water_files/PBS_Files")

script.file<-'/scratch/jc227089/evo-dispersal/art_wbdies/Pilbara/Pilbara_barriers_sims_180days.R'


for (ss in 1:200) {
  
    
    	fid<-ss
	    ##create the sh file
	    zz = file(paste("Pilb_Bar", fid,'.sh',sep=''),'w')
  	  cat('##################################\n',file=zz)
	    cat('#!/bin/sh\n',file=zz)
	    cat('cd $PBS_O_WORKDIR\n',file=zz)
	    cat('source /etc/profile.d/modules.sh\n',file=zz) ### Runs an .sh file which allows modules to be loaded
	    cat('module load R\n', file=zz)
  	  	cat("R CMD BATCH --no-save --no-restore '--args File.ID=", fid, "' ", sep="", file=zz)
	    cat(script.file, " ", paste("Pilb_Bar", fid,'.Rout',sep=''), "\n", sep="",file=zz)
	    cat('##################################\n',file=zz)
  	  close(zz)
			
    	#submit the job
	    system(paste("qsub -m n Pilb_Bar", fid,".sh",sep=""))
  
}

########################
# And then concatenate the output
setwd("~/final_natural_water_files/Pilbara")
point<-c()
times<-c()
nns<-c()
for (ii in 1:200){
  load(file=paste("nn_180", ii, ".RData", sep=""))
  temp<-do.call(what="rbind", output)
  times<-c(times, unlist(temp[,"gen"]))
  point<-c(point, unlist(temp[,"point"]))
  nns<-c(nns, unlist(temp[,"nns"]))
}

summ<-tapply(times<100, list(nns, point), sum)/1000
save(summ, file="Barriers summary 180d.RData")
rm(summ)


pdf(file="n_remove_180Comp.pdf", height=9, width=18)
	par(mfrow=c(1,2))
	load(file="Barriers summary.RData")
	colrs<-c("black", "darkorange", "green")
	nn<-seq(40, 110, 5)

	plot(c(40, 110), c(0, max(summ)), pch=19, xlab="Number of artificial waterbodies removed", 
		ylab="Probability of toads reaching the Pilbara", bty="l", type="n",
		ylim=c(0,1), main="Mean annual rainfall") 
	for (ii in 1:3){
		points(nn, summ[,ii], pch=19, col=colrs[ii])
		lines(nn, summ[,ii], col=colrs[ii])
	}
	legend("topright", legend=c("Eastern", "Central", "Western"), pch=19, col=colrs, bty="n")

	load(file="Barriers summary 180d.RData")
	colrs<-c("black", "darkorange", "green")
	nn<-seq(40, 110, 5)

	plot(c(40, 110), c(0, max(summ)), pch=19, xlab="Number of artificial waterbodies removed", 
		ylab="Probability of toads reaching the Pilbara", bty="l", type="n",
		ylim=c(0,1), main="180 days of movement") 
	for (ii in 1:3){
		points(nn, summ[,ii], pch=19, col=colrs[ii])
		lines(nn, summ[,ii], col=colrs[ii])
	}
	legend("topright", legend=c("Eastern", "Central", "Western"), pch=19, col=colrs, bty="n")

dev.off()






